      subroutine parmovi 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use corgan_com_M 
      use cindex_com_M 
      use numpar_com_M 
      use cprplt_com_M, ONLY:  
      use cophys_com_M 
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx 
      use blcom_com_M, ONLY: drift, qom, t_wall, iplost, plost_pos, plost_neg, &
         iphd2, wate, divpix, bxv, byv, bzv, ex, ey, ez, x, y, z, area_x, &
         area_y, area_z, tsix, tsiy, tsiz, iphead, ijkcell, ijktmp2, ijktmp3, &
         ijktmp4, ijkctmp, nux, nuy, nuz, etax, etay, etaz, px, py, pz, pxi, &
         peta, pzta, up, vp, wp, link, ico, qpar, xptilde, yptilde, zptilde, &
         uptilde, vptilde, wptilde, bxpn, bypn, bzpn, expn, eypn, ezpn, vrms, &
         killer, bcpl, bcpr, bcpt, bcpb, bcpe, bcpf, xl, xr, yt, yb, ze, zf 
      use objects_com_M, ONLY: chi 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv 
      integer , dimension(4) :: lioinput, lioprint 
      integer :: lpx, n, np, is, inner, nout, inew, jnew, knew, ijknew, nsub, isub
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2 
      real(double) :: ixi4, ieta4, ixi5, ieta5, dtosub2 
      real(double), dimension(8) :: wght 
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm 
      real(double), dimension(3,20) :: wsin, wcos 
      real(double), dimension(3,12,20) :: tsi, rot 
      real(double) :: nuxp, nuyp, nuzp, nu, jdotepar, dto2,  rcoll, qomdt&
         , omdtsq, denom, ut, vt, wt, udotb 
      logical :: expnd, nomore, succes 
      character :: name*80 
!-----------------------------------------------
!
!
!
!     *******************************************
!
!     an explicit particle mover for particles in a static
!     magnetic field
!     the equations of motion are solved in natural coordinates
!
!     *********************************************
!
!     jdotepar=0.
!
      lpx = nphist + 1 
      nsub=1
      dto2 = 0.5*dt/dble(nsub)

         do isub= 1, nsub
         
         
!
!dir$ ivdep
!
      iphd2(ijkcell(:ncells)) = 0 
      iphd2(1) = iphead(1) 
!
!
    1 continue 
      nomore = .TRUE. 
!
!     check for more particles to process
!
      ncellsp = 0 
      do n = 1, ncells 
  211    continue 
         if (iphead(ijkcell(n)) <= 0) cycle  
         ijk = ijkcell(n) 
         np = iphead(ijk) 
         is = ico(np) 
         if (.not.drift(is)) then 
!
            ncellsp = ncellsp + 1 
            ijkctmp(ncellsp) = ijkcell(n) 
            ijktmp2(ncellsp) = ijkcell(n) 
!
         else 
!
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijk) 
            iphd2(ijk) = np 
!
            go to 211 
!
         endif 
!
      end do 
!
!     if there are more, process particles as though there were one
!     in every cell in the mesh
!
      if (ncellsp /= 0) then 
!
         nomore = .FALSE. 
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            xptilde(ijk) = px(np) 
            yptilde(ijk) = py(np) 
            zptilde(ijk) = pz(np) 
!
         end do 
!
!
!     **************************************************
!
!
!     corrector step
!
!
!     *************************************************
!
         do inner = 1, 3 
!
            call watev (ncellsp, ijkctmp, iphead, itdim, pxi, peta, pzta, wate) 
!
!
            call trilinp (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, iphead, &
               pxi, peta, pzta, wate, ex, ey, ez, expn, eypn, ezpn) 
!
!			Add constant imposed electric field
!                
            expn=expn+ex_const   
            eypn=eypn+ey_const 
            ezpn=ezpn+ez_const             
!
            call trilinp (ncellsp, ijkctmp, itdim, iwid, jwid, kwid, iphead, &
               pxi, peta, pzta, wate, bxv, byv, bzv, bxpn, bypn, bzpn) 
!
!
!
!     **************************************************************
!
!     solve the momentum equation
!
!     ************************************************************
!
!dir$ ivdep
!
            do n = 1, ncellsp 
!
               ijk = ijkctmp(n) 
               np = iphead(ijk) 
!
!               chip = wate(ijk,1)*chi(ijk+iwid) + (wate(ijk,2)*chi(ijk+iwid+&
!                  jwid)+(wate(ijk,3)*chi(ijk+jwid)+(wate(ijk,4)*chi(ijk)+(wate(&
!                  ijk,5)*chi(ijk+iwid+kwid)+(wate(ijk,6)*chi(ijk+iwid+jwid+kwid&
!                  )+(wate(ijk,7)*chi(ijk+jwid+kwid)+wate(ijk,8)*chi(ijk+kwid)))&
!                  )))) 
!
!               chip = 0. 
!               rcoll = 1./(1. + chip*dto2) 
	rcoll=1d0
!
!
               is = ico(np) 
               qomdt = qom(is)*dto2*rcoll 
!
!
               omdtsq = qomdt**2*(bxpn(ijk)**2+bypn(ijk)**2+bzpn(ijk)**2) 
               denom = 1./(1. + omdtsq) 
!
!
!     solve the momentum equation
!
               ut = up(np)*rcoll + qomdt*expn(ijk) 
               vt = vp(np)*rcoll + qomdt*eypn(ijk) 
               wt = wp(np)*rcoll + qomdt*ezpn(ijk) 
!
               udotb = ut*bxpn(ijk) + vt*bypn(ijk) + wt*bzpn(ijk) 
!
               uptilde(ijk) = (ut + qomdt*(vt*bzpn(ijk)-wt*bypn(ijk)+qomdt*&
                  udotb*bxpn(ijk)))*denom 
               vptilde(ijk) = (vt + qomdt*(wt*bxpn(ijk)-ut*bzpn(ijk)+qomdt*&
                  udotb*bypn(ijk)))*denom 
               wptilde(ijk) = (wt + qomdt*(ut*bypn(ijk)-vt*bxpn(ijk)+qomdt*&
                  udotb*bzpn(ijk)))*denom 
            end do 
!
!dir $ivdep
!
            do n = 1, ncellsp 
!
               ijk = ijkctmp(n) 
               np = iphead(ijk) 
!
               px(np) = xptilde(ijk) + uptilde(ijk)*dto2 
               py(np) = yptilde(ijk) + vptilde(ijk)*dto2 
               pz(np) = zptilde(ijk) + wptilde(ijk)*dto2 
!      px(np)=xptilde(ijk)+uptilde(ijk)*dt*cntr
!      py(np)=yptilde(ijk)+vptilde(ijk)*dt*cntr
!      pz(np)=zptilde(ijk)+wptilde(ijk)*dt*cntr
!     px(np)=xptilde(ijk)
!     py(np)=yptilde(ijk)
!     pz(np)=zptilde(ijk)
!
            end do 
!
            call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
               vrms, t_wall, ico, ijktmp2, ijktmp3, ijktmp4, rmaj, dz, divpix, &
               area_x, area_y, area_z, dto2*2d0, itdim, npart, ibar, jbar, kbar, &
               mgeom, cdlt, sdlt, plost_pos, plost_neg, wate, x, y, z, bxv, byv&
               , bzv, xptilde, yptilde, zptilde, uptilde, vptilde, wptilde, &
               tsix, tsiy, tsiz, etax, etay, etaz, nux, nuy, nuz, link, iplost&
               , qpar, px, py, pz, up, vp, wp, pxi, peta, pzta, cartesian, &
               killer, bcpl, bcpr, bcpb, bcpt, bcpe, bcpf, dx, dy, dz, xl, xr, &
               yb, yt, ze, zf, nu_len, nu_comm) 
 
!	close inner loop
!
         end do 
!
!
!     advance solution to final values
!
!
!dir $ivdep
!
         do n = 1, ncellsp 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
            is = ico(np) 
!
!      jdotepar=jdotepar-.5*qpar(np)*(up(np)**2+vp(np)**2+
!     &   wp(np)**2)/qom(is)
!
            up(np) = 2.*uptilde(ijk) - up(np) 
            vp(np) = 2.*vptilde(ijk) - vp(np) 
            wp(np) = 2.*wptilde(ijk) - wp(np) 
!
!      jdotepar=jdotepar+.5*qpar(np)*(up(np)**2+vp(np)**2+
!     &   wp(np)**2)/qom(is)
!
            px(np) = xptilde(ijk) + uptilde(ijk)*dto2*2d0
            py(np) = yptilde(ijk) + vptilde(ijk)*dto2*2d0
            pz(np) = zptilde(ijk) + wptilde(ijk)*dto2*2d0 
!
         end do 
!
         call parlocat (ncellsp, ijkctmp, iphead, iwid, jwid, kwid, nsampl, &
            vrms, t_wall, ico, ijktmp2, ijktmp3, ijktmp4, rmaj, dz, divpix, &
            area_x, area_y, area_z, dto2*2d0, itdim, npart, ibar, jbar, kbar, mgeom, &
            cdlt, sdlt, plost_pos, plost_neg, wate, x, y, z, bxv, byv, bzv, &
            xptilde, yptilde, zptilde, uptilde, vptilde, wptilde, tsix, tsiy, &
            tsiz, etax, etay, etaz, nux, nuy, nuz, link, iplost, qpar, px, py, &
            pz, up, vp, wp, pxi, peta, pzta, cartesian, killer, bcpl, bcpr, &
            bcpb, bcpt, bcpe, bcpf, dx, dy, dz, xl, xr, yb, yt, ze, zf, nu_len&
            , nu_comm) 
 
!
         nout = ncellsp 
!
!
         do n = 1, ncellsp 
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
!
            inew = int(pxi(np)) 
            jnew = int(peta(np)) 
            knew = int(pzta(np)) 
!
!chek      if((inew-2)*(ibp1-inew).lt.0)  write(*,*) 'np=',np,'inew=',inew
!chek      if((jnew-2)*(jbp1-jnew).lt.0)  write(*,*) 'np=',np,'jnew=',jnew
!chek      if((knew-2)*(kbp1-knew).lt.0)  write(*,*) 'np=',np,'knew=',knew
!
!      brute force
!
            jnew = max(jnew,2) 
            jnew = min(jnew,jbp1) 
            knew = max(knew,2) 
            knew = min(knew,kbp1) 
!
            ijknew = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!     test about locator
!      if(pz(np).lt.z(ijknew).or.pz(np).gt.z(ijknew+kwid)) then
!      write(*,*)'LOCATOR FALLISCE',pz(np),z(ijknew),z(ijknew+kwid)
!      endif



            if (killer(n) /= 0) ijknew = 1 
!
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijknew) 
            iphd2(ijknew) = np 
!
         end do 
!
!     ********************************************************
!
!     a routine to locate a particle on a grid
!     ************************************************************************
!
!     all particles have been processed
!
!     ***********************************************************************
 
         if (.not.nomore) go to 1 
!
      endif 
!
      do n = 1, ncells 
         ijk = ijkcell(n) 
         iphead(ijk) = iphd2(ijk) 
         iphd2(ijk) = 0 
      end do 
      iphead(1) = iphd2(1) 
!
	enddo
!
      return  
!l    ------------------------------------------> return ----->>>
!
      end subroutine parmovi 
